#This script executes event-study regressions to estimate cohort-specific 
#associations between first-generation prenatal exposure to nonattainment 
#and different measures of TSP exposure for first- and second-generation individuals.

rm(list=ls())
gc()

library(reldist)
library(ineq)
library(stargazer)
library(data.table)
library(dplyr)
library(dtplyr)
library(ggplot2)
library(stringdist)
library(lfe)
library(haven)
library(readr)
library(broom)
library(tidylog)
library(rdrobust)
library(HonestDiD)
library(did)
library(here)
library(fixest)

# Set the working directory
setwd("/projects/opp_env/caa1970/")

# Source external R scripts for custom functions
source("Code/Child_Outcomes/child_outcomes_functions.R")
source("Code/rounding.r")

# Load the second-generation analysis dataset
load("Data/second_gen_college_analysis_data_jpemicro.rda")

# Event Studies

#Figure 1a -- First-Generation Prenatal TSP Exposure

event1a <- felm(all_fullterm ~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                
               data=regdata%>% filter(year%in%1969:1975, !is.na(all_fullterm), !is.na(parent_race), !is.na(PI_PC))) 
summary(event1a)

#Save the underlying points
fig1a <- event1a %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(fig1a, file="JPE_Micro/Output/Figure 1/Fig1a.csv", row.names=F)

#Figure 1b - First-Generation Adult TSP Exposure

event1b <- felm(adult_TSP_HS ~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                
               data=regdata%>% filter(year%in%1969:1975, !is.na(all_fullterm), !is.na(parent_race), !is.na(PI_PC))) 
summary(event1b)

#Save the underlying points
fig1b <- event1b %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(fig1b, file="JPE_Micro/Output/Figure 1/Fig1b.csv", row.names=F)

#Figure 1c - Second-Generation Prenatal TSP Exposure

event1c <- felm(child_fullterm ~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                
               data=regdata%>% filter(year%in%1969:1975, !is.na(all_fullterm), !is.na(parent_race), !is.na(PI_PC))) 
summary(event1c)

#Save the underlying points
fig1c <- event1c %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(fig1c, file="JPE_Micro/Output/Figure 1/Fig1c.csv", row.names=F)

#Figure 1d - Second-Generation Pre-College TSP Exposure

event1d <- felm(hs_TSP ~inst_1969+inst_1970+inst_1972+
                 inst_1973+inst_1974+inst_1975 |state_year+county_factor+year_factor+month+
                 survey_year_by_year | 0 |county_factor,                
               data=regdata%>% filter(year%in%1969:1975, !is.na(all_fullterm), !is.na(parent_race), !is.na(PI_PC))) 
summary(event1d)

#Save the underlying points
fig1d <- event1d %>% tidy() %>% mutate_each(funs(signif(.,4)),-term)  %>% 
  filter(grepl("inst", term), !is.na(estimate))%>%
  mutate(year=as.numeric(gsub("inst_","", gsub("_bin", "", term)))) %>% 
  bind_rows(., data.frame(estimate=0, std.error=0, year=1971)) %>% filter(year<=1975) 
write.csv(fig1d, file="JPE_Micro/Output/Figure 1/Fig1d.csv", row.names=F)
